; docformat = 'rst'
;
;+
;
; :Purpose:
;   Aggregate grid cells within local regions.  Aggregated region size depends on the emissions
;     within each region and the average footprint at each grid point.
;     
;   See supplement to ACP paper for discussion on how this is achieved.
;   
; :Inputs:
;   sim_name: (string) Simulation name
;   
; :Requires::
;   regions.sav
;   
;   state.sav
;   
;   LPDM footprints
;   
;   Emissions within each LPDM region
;   
;   If Do_plot is set, David Fanning's Coyote library must be installed
;   
; :Keywords:
;   Do_plot: (input) Set to TRUE to produce encapsulated postscript files showing aggregated areas ({output_directory}/{sim_name}/aggregated_regions_{lri}.eps)
;
; :Outputs:
;   Files defining the grid locations within each aggregated regions
;   ({input_directory}/{sim_name}/ar_{lri}.sav)
;
; :History:
; 	Written by: Matt Rigby, MIT, Aug 19, 2011
;
;-
Pro cels_aggregate, sim_name, do_plot=do_plot

  compile_opt idl2
  on_error, 2

	if keyword_set(do_plot) eq 0 then do_plot=0

  Stations=cels_get_parameter(sim_name, 'STATIONS')
  target=cels_get_parameter(sim_name, 'AGGREGATE_TARGET')
  molmass=cels_get_parameter(sim_name, 'MOLMASS')

  Restore, cels_filestr(/Input, sim_name + '/regions.sav')
  Restore, cels_filestr(/Input, sim_name + '/state.sav')

  nRegions_total=0L
  
  for lri=0, N_LR-1 do begin

    wh_station=where(station_LR eq lri, count)
    if count eq 0 then message, 'No stations in region'

    ;Calculate average footprint
    ;Currently just averages over all files for each region, not necessarily within the date range
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    fname=!null
    for whi=0, n_elements(wh_station)-1 do begin
      fname=[fname, file_search(cels_filestr(/Input, 'fp/'), strcompress('*' + stations[wh_station[whi]] + '*', /remove_all))]    
    endfor
    if fname eq !null then message, strcompress('Footprint files not found for region ' + string(lri, format='(I)'))

    fp_av=fltarr(n_elements(lon_local[lri]), n_elements(lat_local[lri]))
    count=0L
    for fi=0, n_elements(fname)-1 do begin
      restore, fname[fi]
      for ti=0, n_elements(time)-1 do fp_av+=fp[*, *, ti]
      count+=n_elements(time)
    endfor
    fp_av=fp_av/float(count)
    undefine, fp, lon, lat

    ;Calculate average emissions
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    fname=file_search(cels_filestr(/Input, sim_name + '/emissions/'), $
      strcompress('lr_' + string(lri, format='(I02)') + '*', /remove_all)) 
    if fname[0] eq '' then message, strcompress('Emissions files not found for region ' + string(lri, format='(I)'))

    q_av=fltarr(n_elements(lon_local[lri]), n_elements(lat_local[lri]))
    count=0L
    for fi=0, n_elements(fname)-1 do begin
      restore, fname[fi]
      q_av+=q
    endfor
    undefine, q, lonq, latq
    q_av=q_av/float(n_elements(fname))
    q_av=q_av/molmass*1000.*1.e12  ;Convert emissions to pmol/m2/s

    
    ;Define some arrays
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    lon=lon_local[lri]
    lat=lat_local[lri]

    LonSize=n_elements(lon)
    LatSize=n_elements(lat)
    
    lon_2d=fltarr(LonSize, LatSize)
    lat_2d=fltarr(LonSize, LatSize)
    
    for LatI=0, LatSize-1 do Lon_2d[*, LatI]=Lon
    for LonI=0, LonSize-1 do Lat_2d[LonI, *]=Lat
    
  	area=areagrid(lon, lat)

    ;Find 'center' of footprint
  	dummy=max(fp_av, i)
  	ai0=array_indices(fp_av, i)
  	distarr0=mapdist(lon[ai0[0]], lat[ai0[1]], lon, lat, /nowrap)

    ;Find countries in region
    country=find_country(lon, lat, /ocean)
  
    ;Calculate average signal that each element contributes to measurement (ppt)
  	signal=fltarr(LonSize, LatSize)
    signal[lr_wh[lri]]=q_av[lr_wh[lri]]*fp_av[lr_wh[lri]]

    ;Array to contain region numbers
  	region=long(signal gt 0. and finite(country))
  	region[where(region)]=-1L
  	
  	
  	;Begin creating regions
  	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  	ri=1
  	region_total_rec=[0.]
  	count=0L
  	
  	repeat begin
  		region_total=0.
  		region_total+=max(signal*float(region eq -1), i)
  		region[i]=ri
  		count++
      ai=array_indices(signal, i)
      
      dist=distarr0[ai[0], ai[1]]
      
      incr=1
      
  		while (region_total lt target) do begin

        ar_signal=littleregion(signal, ai[0], ai[1], indices=ai_lr, /nowrap, regionsize=max([4, incr^2]))
        ar_country=littleregion(country, ai[0], ai[1], indices=ai_lr, /nowrap, regionsize=max([4, incr^2]))
        ar_region=littleregion(region, ai[0], ai[1], indices=ai_lr, /nowrap, regionsize=max([4, incr^2]))

        wh=where(ar_region eq -1 and ar_country eq country[ai[0], ai[1]])

        if wh[0] ne -1 then begin
          region_total+=total(ar_signal[wh])
          for n=0L, n_elements(wh)-1 do begin
            region[ai_lr[0, wh[n]], ai_lr[1, wh[n]]]=ri
          endfor
          count++
        endif
        
        if (total(region eq -1) eq 0.) then break
        if incr gt round(dist) then break

        incr+=2

  		endwhile
  
      region_total_rec=[region_total_rec, region_total]
  		ri++
  	endrep until (total(region eq -1) eq 0.)
  	ri_max=ri-1

    ;Tidy up small regions
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    wh=where(region_total_rec lt target/10. and region_total_rec gt 0.)
    while wh[0] ne -1 do begin
      region_total=0.
      wh=wh[sort(region_total_rec[wh])]
      ;Find approx center
      wh_region=where(region eq wh[0])
      ai=array_indices(region, wh_region)
      lonI=round(mean(ai[0, *]))
      LatI=round(mean(ai[1, *]))
      dist=mapdist(LonI, LatI, indgen(LonSize), indgen(latsize), /nowrap)
      ri=1.
      repeat begin
        region_ri=region[where(dist lt ri)]
        ri+=0.5
      endrep until total(region_ri ne wh[0] and region_ri ne 0) gt 0.
      region_ri=region_ri[where(region_ri ne 0 and region_ri ne wh[0])]
      region[wh_region]=region_ri[0]
      region_total_rec[region_ri[0]]+=region_total_rec[wh[0]]
      region_total_rec[wh[0]]=0.
      wh=where(region_total_rec lt target/10. and region_total_rec gt 0.)
    endwhile

    
    ;Re-label regions
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    wh=where(region_total_rec gt 0.)
    region_temp=intarr(lonsize, latsize)
    for ri=0, n_elements(wh)-1 do begin
      region_temp[where(region eq wh[ri])]=ri+1
    endfor
    region=region_temp
    ri_max=max(region)

    
    ;Plot regions, if requested
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
  	if do_plot then begin
	  	mr_plot_map, filename=cels_filestr(/Output, sim_name + '/aggregated_regions_' + string(lri, format='(I02)') + '.eps'), $
	  	  float(region), lon, lat, $
        limit=[lat_domain[0, lri], lon_domain[0, lri], lat_domain[1, lri], lon_domain[1, lri]], $
        range=[0., max(float(region))+1.], ct=27, /square, $
        cbtitle='Region index'
		endif

    
    ;Store aggreagted region locations (maybe re-name this to ar_wh?)
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    ar_wh=hash()
    for ri=1, ri_max do begin
      wh=where(region eq ri)
      if wh[0] ne -1 then ar_wh[ri-1]=wh
    endfor

    save, filename=cels_filestr(/Input, sim_name + '/ar_' + string(lri, format='(I02)') + '.sav'), ar_wh
    

    ;Check that emissions still match
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    q_agg=0.
    for n=0, n_elements(ar_wh.Keys())-1 do begin
      if total(q_av[ar_wh[n]]) eq 0. then stop
      q_agg+=total(q_av[ar_wh[n]])
    endfor

    if long(total(q_av[lr_wh[lri]])) ne long(q_agg) then begin
      message, strcompress('Aggregated emissions do not agree agree: ' + $
        string(total(q_av[lr_wh[lri]])) + ' ' + string(q_agg))
    endif
    
    nRegions_total+=n_elements(ar_wh->keys())

  endfor

  print, 'CELS_AGGREGATE: Total aggregated regions ', nRegions_total

End